<!DOCTYPE html>
<html>
<head>
<link href='https://fonts.googleapis.com/css?family=Roboto Condensed' rel='stylesheet'>
<meta name="viewport" content="width=device-width, initial-scale=1.0">

<style>

body {
	margin-top: 1cm;
	margin-left: 4cm;
	margin-right: 4cm;
	padding: 0;
        font-size: 15px;
        font-family: 'Roboto Condensed', 'Tahoma', 'Arial', sans-serif;
}

   ul.tree a {
    color: back;
    text-decoration: none;
    display: block;
    width: 100%;
   }
   ul.tree a:hover {
    background: lightgray;
    text-decoration: none;
    display: inline-block;
    width: 100%;
   }

   li {
     list-style-type: none;
     background:  repeat-y;
     padding: 0;
   }


   ul.tree, ul.tree ul {
     list-style-type: none;
     background:  repeat-y;
     margin: 0;
     padding: 0;
     margin: 5px;
   }

   ul.tree ul {
     list-style-type: none;
     margin-left: 10px;
   }


a {
    color: black;
}

   ul.tree li {
     margin: 5px;
     padding: 0 12px;
     line-height: 20px;
     color: #369;
   }

   ul.tree li.last {
     background: #fff url(images/lastnode.png) no-repeat;
   }

ul.topnav {
    list-style-type: none;
    padding: 0;
    margin: 0;
    overflow: hidden;
    background-color: #333;
    font-size: 15px;
    text-align: center;
}

ul.topnav li {float: left;}

img.logo {max-width: 350px;}

ul.topnav li a {
    display: block;
    color: white;
    text-align: center;
    padding: 14px 27px;
    text-decoration: none;
}
ul.topnav li a:hover:not(.active) {background-color: #111;}
ul.topnav li a.active {background-color: #4CAF50;}
ul.topnav li.right {float: right;}
@media screen and (max-width: 600px){
    ul.topnav li.right, 
    ul.topnav li {float: none;}
}

pre, code { font-size:12px; }
tt { font-size: 1.2em; }
pre { margin:0px 0px 20px; }
pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }
pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }
pre.error { color:red; }

@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }

span.keyword { color:#0000FF }
span.comment { color:#228B22 }
span.string { color:#A020F0 }
span.untermstring { color:#B20000 }
span.syscmd { color:#B28C00 }

.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }
.footer p { margin:0px; }
.footer a { color:#878787; }
.footer a:hover { color:#878787; text-decoration:underline; }
.footer a:visited { color:#878787; }

table th { padding:7px 5px; text-align:left; vertical-align:middle; border: 1px solid #d6d4d4; font-weight:bold; }
table td { padding:7px 5px; text-align:left; vertical-align:top; border:1px solid #d6d4d4; }

</style>
<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_CHTML"></script>
</head>
<body> 
<img class="logo" src="ValoMC_logo.png">
<ul class="topnav">
<font size="4">
  <li><a href="index.html">Home</a></li>
  <li><a href="download.html">Download</a></li>
  <li><a href="features.html">Features</a></li>
  <li><a href="installation.html">Installation</a></li>
  <li><a class="active" href="documentation.html">Documentation</a></li>
  <li class="right"><a href="about.html">About</a></li>
</font>
</ul>

<div style="padding:0 16px;">
  <h2></h2>
  <p></p>
  <p></p>
<font size="4">
<div style="width: 100%;">
    <div style="width: 300px; float: left;">
<br> 
<li><a>Documentation</a></li>
<ul class="tree">
<li><a href="gettingstarted.html">Getting started</a></li>
<li><a>2D</a></li>
<ul class="tree">
<li><a href="simpletest.html">Simple example</a></li>
<li><a href="inhomogeneous.html">Creating an inhomogeneous medium</a></li>
<li><a href="directingls.html">Directing light sources</a></li>
<li><a href="boundarytest.html">Setting boundary conditions and visualising the boundary solution</a></li>
<li><a href="frequency.html">Frequency domain calculation</a></li>
<li><a href="toasttest.html">Working with Toast</a></li>
<li><a href="netgentest.html">Working with NetGen</a></li>
<li><a href="generatingc.html">Generating input for the external executable</a></li>
<li><a href="pixeltest.html">Working with pixel format data</a></li>
<li><a href="kwavetest.html">Simulating the photoacoustic effect using K-Wave</a></li>
<li><a>Code documentation</a></li>
<ul class="tree">
<li><a href="structures.html">List of structures</a></li>
<li><a href="findingelements.html">Finding elements</a></li>
<li><a href="findingboundaries.html">Finding boundaries</a></li>
</ul></ul>
<li><a>3D</a></li>
<ul class="tree">
<li><a href="voxeltest.html">Working with voxel format data</a></li>
<li><a href="netgentest3d.html">Working with NetGen</a></li>
<li><a href="threedmodel.html">Digimouse example</a></li>
<li><a href="kwavetest3d.html">Simulating the photoacoustic effect using K-Wave</a></li>
<li><a>Code documentation</a></li>
<ul class="tree">
<li><a href="structures3d.html">List of structures</a></li>
<li><a href="findingelements3d.html">Finding elements</a></li>
<li><a href="findingboundaries3d.html">Finding boundaries</a></li>
</ul></ul>
<li><a href="literature.html">Literature</a></li>
<li><a href="functionreference.html">Alphabetical function listing</a></li>
</ul>
    </div>
    <div style="margin-left: 350px; padding: 5px;"> 
<div class="content"><h2>Contents</h2><div><ul><li><a href="#1">Simulating the photoacoustic effect using K-Wave</a></li><li><a href="#2">Create the k-Wave grid</a></li><li><a href="#3">Create a PhotonMC mesh</a></li><li><a href="#4">Define optical coefficients</a></li><li><a href="#5">Create a light source</a></li><li><a href="#6">Compute the initial pressure from the photon fluence</a></li><li><a href="#7">Define the k-Wave sensor mask</a></li><li><a href="#8">Move the perfectly matched layer (PML) outside of the computation domain and run the acoustic simulation</a></li></ul></div><h2 id="1">Simulating the photoacoustic effect using K-Wave</h2><p>This example demonstrates simulation of a pressure field generated through the absorption of an externally introduced light pulse. The light propagation is simulated using PhotonMC and the propagation and detection of pressure wavefield is simulated using k-Wave, see <a href="http://www.k-wave.org/documentation/k-wave_initial_value_problems.php">http://www.k-wave.org/documentation/k-wave_initial_value_problems.php</a>. The example also shows how the computation grid of k-Wave and mesh of PhotonMC can be made compatible. Note that k-Wave uses SI units (e.g. [m]) and PhotonMC works in millimetre-scale (e.g. [mm]). k-Wave is an open source acoustics toolbox for MATLAB and C++ developed by Bradley Treeby and Ben Cox (University College London) and Jiri Jaros (Brno University of Technology). The software is designed for time domain acoustic and ultrasound simulations in complex and tissue-realistic media.</p><p>k-Wave homepage: <a href="http://www.k-wave.org/">http://www.k-wave.org/</a>. B. E. Treeby and B. T. Cox: "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields", Journal of Biomedical Optics, 15(2):021314, 2010.</p><pre class="codeinput">clear <span class="string">all</span>;
</pre><h2 id="2">Create the k-Wave grid</h2><pre class="codeinput">Nx = 60;           <span class="comment">% number of grid points in the x (row) direction</span>
Ny = 60;           <span class="comment">% number of grid points in the y (column) direction</span>
Nz = 60;           <span class="comment">% number of grid points in the y (column) direction</span>
dx = 0.1e-3;        <span class="comment">% grid point spacing in the x direction [m]</span>
dy = 0.1e-3;        <span class="comment">% grid point spacing in the y direction [m]</span>
dz = 0.1e-3;        <span class="comment">% grid point spacing in the y direction [m]</span>
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

<span class="comment">% Create two internal structures using makeDisk</span>
sphere = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, 5);


<span class="comment">% Define the acoustic properties</span>
sphere_indices = find(sphere==1);

medium.sound_speed = 1500*ones(Nx, Ny, Nz);    <span class="comment">% [m/s]</span>
medium.sound_speed(sphere_indices) = 1800;     <span class="comment">% [m/s]</span>

medium.density = 1000*ones(Nx, Ny, Nz);        <span class="comment">% [kg/m^3]</span>
</pre><pre class="codeoutput">WARNING: makeGrid will be deprecated in a future version of k-Wave.
         Update codes to use the syntax kgrid = kWaveGrid(...).
</pre><h2 id="3">Create a PhotonMC mesh</h2><p>PhotonMC uses triangles and tetrahedrons as the basis shape, whereas in k-Wave pixels and voxels are used. createGridMesh can be used to create a straightforward mapping between the triangles and pixels. Note that K-Wave uses matrices in the format matrix(X,Y,Z), whereas MATLAB t(c.f. meshgrid) and PhotonMC uses matrix(Y,X,Z) Therefore x and y should be swapped when moving between PhotonMC arrays and k-Wave arrays.</p><pre class="codeinput">pmcmesh = createGridMesh(kgrid.y_vec*1e3, kgrid.x_vec*1e3, <span class="keyword">...</span>
                         kgrid.z_vec*1e3); <span class="comment">% [m to mm]</span>
</pre><h2 id="4">Define optical coefficients</h2><p>For users accustomed to k-Wave, the optical coefficients can be set in similar fashion as in k-Wave, i.e. using multidimensional arrays. If one of the arrays defining the medium is given as multidimensional array to PhotonMC, the code will assume that the mesh was created using 'createGridMesh' and the output fluence will also given as a two dimensional array in solution.grid_fluence.  See the example 'Working with pixel and voxel data' on how to achieve similar assignments using one dimensional indexing.</p><pre class="codeinput">pmcmedium.scattering_coefficient = 0.01*ones(Nx, Ny, Nz);
pmcmedium.absorption_coefficient = 0.001*ones(Nx, Ny, Nz);

pmcmedium.absorption_coefficient(sphere_indices) = 0.02;
pmcmedium.scattering_anisotropy = 0.9;           <span class="comment">% scattering</span>
                                                 <span class="comment">% anisotropy</span>
                                                 <span class="comment">% parameter</span>
                                                 <span class="comment">% [unitless]</span>
pmcmedium.refractive_index = 1.0*ones(Nx, Ny, Nz);

<span class="comment">% Define the Gruneisen parameter describing photoacoustic</span>
<span class="comment">% efficiency</span>
pmcmedium.gruneisen_parameter = 0.02*ones(Nx, Ny, Nz);      <span class="comment">% [unitless]</span>
</pre><h2 id="5">Create a light source</h2><p>Set a light source with a width of 2 mm and cosinic directional profile at the bottom of the computation domain</p><pre class="codeinput">boundary_with_lightsource = findBoundaries(pmcmesh, <span class="string">'direction'</span>, <span class="keyword">...</span>
                                           [0 0 0], <span class="keyword">...</span>
                                           [-5 0 0], <span class="keyword">...</span>
                                           1);

pmcboundary.lightsource(boundary_with_lightsource) = {<span class="string">'cosinic'</span>};

<span class="comment">%%Run the Monte Carlo simulation</span>

solution = ValoMC(pmcmesh, pmcmedium, pmcboundary);
</pre><pre class="codeoutput">Initializing MC3D...
Computation uses 16 threads
Computing... 
...done

Done
</pre><h2 id="6">Compute the initial pressure from the photon fluence</h2><p>Note that since the medium was defined as three dimensional arrays, the output is also given as a three dimensional array.</p><pre class="codeinput"><span class="comment">% Compute the absorbed optical energy density</span>
<span class="comment">% 1e3 converts [J/mm^2] to [J/m^2]</span>
pmcmedium.absorbed_energy = pmcmedium.absorption_coefficient .* <span class="keyword">...</span>
    solution.grid_fluence*1e3; <span class="comment">% [J/m3]</span>

<span class="comment">% Compute the initial pressure distribution</span>
source.p0 = pmcmedium.gruneisen_parameter .* pmcmedium.absorbed_energy; <span class="comment">% [Pa]</span>
</pre><h2 id="7">Define the k-Wave sensor mask</h2><pre class="codeinput"><span class="comment">% Define a circular sensor</span>
sensor_radius = 2e-3;         <span class="comment">% [m]</span>
num_sensor_points = 1000;     <span class="comment">% number of sensor points</span>

sensor.mask = makeCartSphere(sensor_radius, num_sensor_points);
</pre><h2 id="8">Move the perfectly matched layer (PML) outside of the computation domain and run the acoustic simulation</h2><p>The PML is a layer that absorbs waves for simulating free regions and is normally contained within the computation region of k-Wave. For a more straightforward mapping between the computation region of k-Wave and PhotonMC, the PML is moved outside of the computation region.</p><pre class="codeinput">sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, <span class="keyword">...</span>
                                <span class="string">'PMLInside'</span>, false);


[X Y Z] = meshgrid(kgrid.x_vec*1e3, kgrid.y_vec*1e3, kgrid.z_vec* <span class="keyword">...</span>
                   1e3);
slice(X, Y, Z, source.p0, 0, 0, 0);
xlabel(<span class="string">'y [mm]'</span>);
ylabel(<span class="string">'x [mm]'</span>);
zlabel(<span class="string">'z [mm]'</span>);
view(-33,35);
title(<span class="string">'Initial pressure'</span>);

<span class="comment">% plot the simulated sensor data</span>
figure;
imagesc(sensor_data,  [min(sensor_data(:)) max(sensor_data(:))]);
colormap(getColorMap);
ylabel(<span class="string">'Sensor Position'</span>);
xlabel(<span class="string">'Time Step'</span>);
c = colorbar;  <span class="comment">% create a colorbar</span>
colorbar;
title(<span class="string">'Sensor data'</span>);
</pre><pre class="codeoutput">Running k-Wave simulation...
  start time: 23-May-2018 18:38:43
  reference sound speed: 1800m/s
  dt: 16.6667ns, t_end: 6.9167us, time steps: 416
  input grid size: 60 by 60 by 60 grid points (6 by 6 by 6mm)
  maximum supported frequency: 7.5MHz
  smoothing p0 distribution...
  expanding computational grid...
  computational grid size: 80 by 80 by 80 grid points
  calculating Delaunay triangulation...
  precomputation completed in 29.4521s
  starting time loop...
  estimated simulation time 21.0806s...
  simulation completed in 20.5736s
  total computation time 50.1759s
</pre><img alt="" hspace="5" src="kwavetest3d_01.png" vspace="5"/> <img alt="" hspace="5" src="kwavetest3d_02.png" vspace="5"/> <p class="footer"><br/><a href="http://www.mathworks.com/products/matlab/">Published with MATLAB® R2016b</a><br/></p></div>
    </div>
</div>
</font>

  </div>

<br>
<br> </div><footer> <hr> <p>Last updated: Mon May 28 11:28:14 EEST 2018 by aleksle</p></footer> </body> </html>

